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ABSTRACT 

We use the Made-to-Measure method to construct iV-body realizations of self-similar, 
triaxial ellipsoidal halos having cosnrologically realistic density profiles. Our imple¬ 
mentation parallels previous work with a few numerical refinements, but we show that 
orbital averaging is an intrinsic feature of the force of change equation and argue that 
additional averaging or smoothing schemes are redundant. We present models hav¬ 
ing the Einasto radial mass profile that range from prolate to strongly triaxial. We 
use a least-squares polynomial fit to the expansion coefficients to obtain an analytical 
representation of the particle density from which we derive density contours and ec¬ 
centricity profiles more efficiently than by the usual particle smoothing techniques. We 
show that our iV-body realizations both retain their shape in unconstrained evolution 
and recover it after large amplitude perturbations. 

Key words: methods: N-body simulations — methods: numerical — galaxies: halos 
— galaxies: kinematics and dynamics — galaxies: structure 


1 INTRODUCTION 


Dark-matter halos that form hierarchically in simulations 
of ACDM cosmology have a number of regular properties: 
Their mass density profiles, the distributions of velocity 
anisotropy, spin parameters, and substructure are largely 
independent of the halo’s mass over a wide range of val¬ 
ues. Other properties, such as the concentrat ion an d tri- 
axiality, vary systematically with halo mass l|Zem 3 Hoi; 
iDiemand fc MoorelbOlltTFrenk fc W hitell20 12l). 

The landmark work of iNavarro et alJ ( 19961 . hereafter 
NFW) reported that the halo mass density could be fitted 
with a broken power-law 




7-3 


(i) 


where the length scale, r s , is the radius around which the 
logarithmic slope changes. The steepness of the inner cusp 
reflected in the va lue of the exponent 7 has been t he su bject 


of much debate (iMoore et al. 19991: Fukushi 

le fc Makinol 

200ll; iPower et al. 2003[; Diemand et al. 

2004 

20051). More 

recent work (Navarro et al. 2004|, 201C)|: 

Merritt et al.l 20051: 


iGao et al.ll2008l 'l finds the Einasto profile (described in §4) 
with its steep but finite central density to be a better fit. 


Virialized halos that have undergone multiple mergers 
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acquire a triaxial shape i n general, though with a strong 
bias t owards prolateness jFrenk et al.l 1 19881 : ICole fc Lacevl 
1 19961 : iBailin fc Steinmetd l2005li . To accommodate tri axial- 
ity when fitting the mass density, a number of authors l|Kat d 
Il99ll : iDubinski fc Carlberd Il99ll : Ijing fc Sutol 120021 . here¬ 
after JS02) replace the radial coordinate r with the dimen¬ 
sionless coordinate £ defined through 


f 




a ^ b ^ c 


( 2 ) 


where a, b, and c are the lengths of the major, in¬ 
termediate, and minor axes, respectively. Although ax¬ 
ial ratios and alignments di s play some radial varia¬ 
tion (JS 02 : ICole fc Lacevl Il996l : iBailin fc Steinmetd 120051 : 
ISchneider et al. 20121 . h^eafter SFC 12 ), we focus here on 
creating self-similar ellipsoids for which these values remain 
in a fixed ratio, with the axes aligned, at all radii. This ide¬ 
alized geometry reduces the mass density to a function of a 
single variable £ (eq. O, with the axis ratios b/a and c/a as 
constant parameters. 

Self-similarity is an attractive assumption, but it is 
generally insufficient to enable a derivation of a distribu¬ 
tion function ( DF). The DF can be exp ressed in terms 
of the actions teinnev fc Tremaine! boost ). but these are 
known only in the case of the “perfect” ellipsoid for which 
th e potential is fully separable in ellipsoidal coordinates 
(Ide Zeeuw fc Lvnden-Bellll 19851 : Ide Zeeuwlll985l 'l and all or¬ 
bits are regular. Regular orbits also exist in more general 
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2 J. C. Malvido & J. A. Sellwood 


ellipsoids but irregular orbits, for which the only isolat¬ 
ing integral is the energy, popul ate some parts of phase 
space, and may even dominat e dUdr v fc Pfennigen 1 19881 : 
IValluri fe Merril~i1ll998l : IValluri et alj|2010h . 

Various approaches have been adopted to construct N- 
body models of triaxial mass d istributions with some suc - 
cess. iHolIev-Bockelmann et ahl ( 2001) and Widrowl (2008) 

[ 200 a. 


“squeeze” a spherical model, ' [Rodionov et al.l (|20 oA 120091 1 


developed an iterative method that corrects the density after 
each time step, without revising the velocit ies, until equlib- 
rium is attained, while iMoore et al.l (}2004j ) construct them 
through multiple collisions of separate systems. However, 
neither the squeeze nor collision methods provide a system¬ 
atic means for reproducing a specific geometry, such as self¬ 
similarity, while only halos with a central core, rather than a 
steep inner density gradient, have been attempted with the 

iter ative method. _ _ 

ISchwarzschildl d 19791 . 1 19931 ) built triaxial models from 
linear combinations of bound orbits in the desired po¬ 
tential. In his method, orbits are drawn from a pre¬ 
calculated library that approximately spans the phase- 
space of the bound system and contains the time-averaged 
contributions to the mass density in a set of spatial 
cells. The assignment of orbital weights is constrained 
by the targeted mass density in these cells, but since 
the assignment is far from unique, it is customary to 
seek an optimal solution, such as by maximizing the en¬ 
tropy (iRichstone fc Tremaine 119881 1 to favor a smooth dis- 


tribution. The method has been successfully apr 

rlied to 

cusdv, triaxial svstems (Merritt & Fridman 19961: 

Merritt] 

19971: Caouzzo-Dolcetta et al. 

200?!: van de Yen et al. 20081: 

Vasiliev & Athanassoula 2012 

). However, creation 

of N- 


body models requires the further, non-trivial step of select¬ 
ing particles from appropriately chosen points along the or¬ 
bits selected from the library. 

Here we adopt the “ made to measure” (her eafter 
M2M) algorithm proposed bv ISver fe Trema ine (1996, here¬ 
after ST96), which is similar in spirit to Schwarzschild’s 
but adjusts the mass of each moving particle dynam- 
ically. The algo r ithm has been further develop e d by 


Dc Lorenzi et al 


hereafter D09), 


ltnm has been lurther develop e d by 
(|2007|, here after D L07), Dehneg_|2003, 
Long fe Mad d201Ch . iMorganti fe Gerhard 


< 20121 ) . and iHunt fe Kawat al <201 fil l. and a useful brief re¬ 
view is given by iGerhardl (1201 (l b A further advantage of 
M2M is that, on completion, the equilibrium V-body sys¬ 
tem is ready to use with no additional step. The range of ap- 
plications of M2M ( |De Lorenzi et ahl 20081 . 120091 ; iDas et ahl 
l201ll ; iLong fe Mad 1201 J Long et aTll2013li is a fair measure 
of its versatility. 

One can add observational constraint s to bo th M2M 
6DL07) and Schwarzschild’s method (see e.g. lChaname et alJ 
120081 . and references therein) in order to construct models 
for comparison with data, which is probably the most widely 
used application of these methods. But our objective in this 
paper is purely theoretical. 

Our aim is to use M2M to construct A-body realizations 
of self-similar ellipsoidal halos with a mass density and ec¬ 
centricity suggested by cosmological simulations. While an 
improvement on spherical halos, our models are still idealiza¬ 
tions, since cosmological halos are generally not self-si milar, 
and we also neglect rota tion fe.g. iBu llock et al .Il 200 l l). cer- 
tain kinematic features (I Cole fe Lacey 19961 : Woitak et al.l 


2005 


2005 


Hansenl 2009lj and s ubstructure llDiemand et al .1120041 . 


Springel et ajj f2008). 


Some triaxial models have already been created by this 
technique, notably in DL07 and D09, but there are minor 
differences in our approach and we present more extensive 
tests of the stability of the resulting models. We will use 
them in later work to study the influence of halo triaxiality 
on the dynamics of embedded disks. 


2 THE ALGORITHM 

The M2M method has already been described elsewhere 
(e.g., ST96, DL07, D09). We therefore give only an outline, 
noting where our treatment differs from previous work. Our 
presentation follows closely that given by DL07, but we do 
not require their generalization to include constraints from 
velocity data. 


2.1 Preliminaries 


We seek to create equilibrium halo models with a target den¬ 
sity distribution pr(x) and total mass My in the gravita¬ 
tional potential $t(x). In most applications, the potential 
arises from the density pr only, but it is also possible to 
consider models with additional mass components, such as 
a disk or central mass, that contribute to but not to pT- 
First we divide the volume of the target model into a 
set of K cells, each having a volume 14, and define a set of 
occupancy kernels {ICk(x), k = 1,..., K} 


f 1 a: is within 14 
L 0 otherwise. 


( 3 ) 


The mass in the fcth cell is therefore given by 


Pk= K.k{x)p T (x)d 3 x. 


( 4 ) 


Suppose that we were to construct a trial set of particles 
(e.g., as described in %L7l) 

N 

To(z) = m p ^^w 0 iS(z -Zi(t = 0)), (5) 

i= 1 

where z = (x,v) is a 6 D phase space coordinate. The initial 
density is determined by the adopted set of N particle posi¬ 
tion coordinates (0)} with their initial fractional weights 
{rooi} normalized by the fiducial particle mass m p (see *13.411 
such that m P Y2 = i w oi = Mt■ At t = 0, this initial guess 
therefore has the mass in each of K cells 

N 

P 0 k = m p ^ woiJCk {Xi (0)). ( 6 ) 

i= 1 

The density at later times, p(x, t), is determined by integrat¬ 
ing all particle trajectories to time t in the frozen potential 
<I>t, when the mass in each cell becomes 

N 

p k (t) = m p ^ Wi(t)JC k (xi(t)). ( 7 ) 

i= 1 

Schwarzschild’s approach is first to compute the time 
average occupancy of each particle in each cell 
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Self-similar triaxial halos 3 


K-ik = - / lCk(Xi(t)) dt, (8) 

T Jo 

for some large but finite time period t, and then to require 
that 

N 

(p k ) = m p ^ WiJCik = Pk , (9) 

i= 1 


for a set of time-independent and non-negative {wi}. Since 
the space of possible solutions to eq. © is generally large, 
it is desirable to optimize for smoothess; a popular choice is 
to maximize the “entropy”, which is defined as 


5 



i=l 


( 10 ) 


We use our initially guessed weights for the priors, i.e. un = 
wot, but other choices are possible. An entropy decrease, 
A<S < 0, may be interpreted as an increase of “order” in the 
physical system; for example, a sphere is less ordered than 
an ellipsoid. 


2.2 Made to measure 


In the M2M method, the particle weights Wi{t) are adjusted, 
while maintaining the condition p(x,t) « pr(x), as the par¬ 
ticle trajectories are integrated. It is usually cast as an op¬ 
timization problem in which the “likelihood function” 

£(f) = pS(t) — C(t), with p > 0 (11) 


is to be maximized. Both parts of this objective are time- 
dependent, since the weights are adjusted. The “cost func¬ 
tion” is 


C(t) 



k = 1 


Pk it) - P k 

0~k 


( 12 ) 


where a k is the dispersion of pok, i.e., the sampling uncer¬ 
tainty implicit in generating Po(z) (DL07, D09). 

The c ondition for the optimal solution for {w *} when 
w* > 0 is dAvriell[l976l i 

dc 

Wi -— = 0. i = 1,..., N (13) 

OWi 

Therefore, either w* = 0, or 


dC dS 1 A K k {Xi) A 
dwi ^ dwi N /j k k 


= 0 , 


(14) 


where A k = ( p k - P k )/a k - 


2.3 The FOCE 

Originally proposed by ST96, the force-of-change-equation 
(FOCE) 

dC 

Wi=eWi -—, (15) 

OWi 

with e being a small parameter, provides a rule for adjusting 
Wi in order to drive dC/dwt towards zero, a procedure in a 
simulation that is readily combined with the particle motion. 

The Wi(t ) and A k (t) functions have noisy paths owing 
to IV-body sampling noise and the movement of particles 


across cells, which causes abrupt changes to the values of 
1Ck(Xi). ST96 suggested replacing the A k in eq. (fl5l) with its 
moving average, Ak, such that A*,(0) = A k (0) and 

^= v (A k -Ak) (16) 

where p ^ 2 e, as a means of suppressing large fluctuations 
in the FOCE from these sources of noise and likened its 
effect to that of smearing a set of virtual particles over a 
given orbit. This adjustment has been integral to imple¬ 
mentations of M2M reported in the literature. D09 points 
out, however, that this procedure is inconsistent with the 
extremum condition m and compromises the meaning of 
C in eq. m since the &k reflect the original A-body shot 
noise, not the reduced noise. Instead, D09 applies a mov¬ 
ing average to dC/dwi, and, in combination with eq. (1151) . 
obtains a second order ODE for In Wi. 

However, the following argument suggests that such at¬ 
tempts at explicit averaging are redundant. Formally, the 
FOCE may be integrated to give 

Wi{t) = e eDi(i) Wio, (17) 

where 

Di(t)= [ dt' . (18) 

Jo 9wi 

Thus we see that eDi ( t ) is proportional to the time average 
{•dC/dwi)t ■ Substitution of eq. (fill) for the integrand in (TTH1) 
shows that the behavior of eDi(t) is dominated by orbital 
averages of the form ( ICk(Xi)A k ) t . But K-k(Xi ) is influenced 
by the trajectory of the ith particle only, while A*, is de¬ 
termined by all particles in the fcth cell. If the number of 
such particles is sufficiently large, it is reasonable to expect 
that ( K,k{Xi)Ak) t decouples at large times into the product 
of two independent averages ICik x (A k ) t , where )Cik is de¬ 
fined in eq. ©• Furthermore, ST96 and DL07 prove that 
(A k) t —> 0 at large times which, from the definition of A*,, 
implies (pk)t —¥ Pk- Thus we effectively recover eq. © and 
thereby formally establish the asymptotic equivalence of the 
M2M and Schwarzschild’s methods. 

Therefore, not only is explicit averaging redundant but 
the application of eq. uni is counterproductive if N is large 
enough. We have not explored how large N needs to be 
for this statement to be true, but when we employ N = 
1.2 x 10 6 particles, we find that inclusion of the averaging 
procedure m degrades the actual value of C by as much as 
an order of magnitude. We therefore refrain from including 
any additional averaging. 

Thus we find a deeper resemblance between the M2M 
approach and Schwarzschild’s method. However, there are 
two reasons for preferring M2M: not only does eq. com¬ 
bine averaging and solving for {w*} into a single algorithm 
instead of two distinct steps but, upon convergence, it also 
delivers a representative set of particle coordinates that are 
in equilibrium, in readiness for an A-body simulation. 


2.4 Convergence 

Fluctuations in the values of the Wi(t) as the system evolves 
jointly with the FOCE imply that satisfying condition ©3 
for all N particles simultaneously at any given instance is 
implausible. Instead, an equally effective, yet more realistic, 
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goal of the M2M method is to ensure that max ||(Afc)[| < 5, 
for some prescribed threshold <5. 

By linearizing the FOCE (exclusive of the entropy 
term), ST96 (see also DL07) obtained a set of equations for 
the modes of the (A*,) and demonstrated that these modes 
decay exponentially - a result that essentially formalizes the 
intuition behind the intent of the FOCE. Though it does 
not provide practical guidelines for assessing convergence, it 
does, however, render the (A*,), rather than the Wi(t), as 
the focus of the convergence analysis. 

The threshold, 5, may be made as small as desired by 
reducing /r, with the disadvantage that an increasing num¬ 
ber of particle weights are driven toward zero. Reaching a 
desired threshold <5 is not necessarily indicative of conver¬ 
gence; only when changes in the entropy become small rel¬ 
ative to its level can the algorithm be said to have reached 
convergent conditions. 

Ultimately, the true test of the algorithm’s convergence 
is whether the resulting system maintains the desired shape 
and mass profile as it evolves freely and self-consistently 
without intervention from the FOCE (see < 15.21) . 

2.5 Summary of the Procedure 

Our implementation of the algorithm, which will be dis¬ 
cussed in greater detail later, may be conveniently summa¬ 
rized as follows (see also D09): 

(i) Initialization : Construct an equibrium, spherically 
symmetric model with the desired mass density profile con¬ 
sisting of Nm particles. Transform the particle coordinates 
by rescaling (see eq. [ 21 ] below) to achieve the targeted ge¬ 
ometry and adjust the velocities so that the tensor virial 
theorem is satisfied. 

(ii) Relaxation: Allow the system to settle (without the 
FOCE) for 20 - 25 dynamical times. 

(iii) M2M step : Evolve the system while adjusting the 
weights in accordance to the FOCE. 

(iv) Convergence check: Once the entropy level is nearly 
constant, switch off the FOCE and allow the system to 
evolve self-consistently to test whether it is in equilibrium. 

(v) Stability Check: Apply an adiabatic, aspherical per¬ 
turbation to seed any possible secular instabilities and follow 
the subsequent self-consistent evolution. 

Note that we obtain a smooth target potential and the 
values of dk needed for steps (ii) and (iii) by creating a model 
having many more particles than is practical to evolve, i.e., 
Nt 3> Nm- Their positions, with velocities ignored, give a 
smoother representation of pr, and we can solve for their 
gravitational field using our adopted A-body method. 

It is possible to run M2M using a self-consistent, rather 
than a rigid, potential. This is less satisfactory, however, 
because M2M continuously tries to adjust weights to com¬ 
pensate for (mild) particle scattering. We therefore prefer to 
use a smooth and frozen target potential. 


3 COMPUTATIONAL METHODS 
3.1 Self-Similar Ellipsoid 

The target triaxial halos in this work consist of a series of 
concentric and co-eccentric, nested ellipsoidal shells - a sim¬ 


plifying idealization. Surfaces of constant density in a self¬ 
similar ellipsiod have constant £, where 


2 2 
>2 2 2 ® 2 ® 

€ =* +y~ ¥ + z 


(19) 


Here a, b, and c are the semi-axes of the outer surface 
of the finite ellipsoid. We also define the eccentricities, 
Sy — (1 — b 2 /a 2 ) 1 / 2 and e z = (1 — (? /a 2 ) 1 ^ 2 for later use. 
Strongly triaxial systems have values of the triaxiality pa¬ 
rameter, T = (a 2 — b 2 )/(a 2 — c 2 ) ~ 0.5, while T = 0 for an 
oblate, and T = 1 for a prolate, spheroid. 

We wish to construct self-similar ellipsiods for which the 
density profile p e (£) has the same functional form (to within 
a normalization factor to preserve the same total mass) as 
some spherical model p s (r). We relate a point r = (x, y , z) on 
a sphere to r' = (x' , y' , z') on the ellipsoidal surface through 


2 2 , 2 2 /2 12 2 /(2 
r =x +y +z =x +y a /b 


/2 2 
+ z a 


a 2 =e 


( 20 ) 


which is equivalent to the scale transformation 
x = x ; y = ya/b ; z = za/c. ( 21 ) 


Two ellipsoidal surfaces at £ and £ + <5£ bound a volume 
known as a thin homoeoid (BT08, §2.5). The volume of a 
thin homoeoid is 47 r d yz ^ 2 d^, with 9 yz = bc/a 2 , while that of 
a spherical shell is 4ivr 2 dr. Therefore, if a sphere of radius a 
and density p s {r) is compressed along the intermediate and 
minor axes by factors b/a and c/a respectively in such a way 
that equidensity shells do not cross, it yields a self-similar 
ellipsoid of equal mass with density p e {x) = 9~f ps{(/)- 


3.2 Force computation 


Particle-mesh, or grid, methods remain the most efficient 
for A-b od y sim ulations of an isolated gravitating system 
(ISellwoodll2014l l. Here we use a radial grid that computes 
forces from an expansion in surface harmonics on a set of 
spherical shells that are more closely spaced near the center. 
An arbitrary field function 4/(x) (e.g., the gravitational po¬ 
tential, force components, mass density, etc.) can be written 
in real form as 

OO l 

*(*) = 

1=0 m=0 

i m {r) cos mcj)+'I’t-mir) sin m^], ( 22 ) 
where x = ( r , 9, cj>) is the field point, 


n ”( s > s < 23 > 

P; m (cosf?) is the associated Legendre function, and the nor¬ 
malization factor 7 i m = (2 — < 5 m o)( 2 Z+l)/ 47 r. The coefficients 
in eq. ( 1221 ) can be computed from 





f cos m</> 
l sin mf> 


4t(r, O ', <j>) sin 9' d9'd(j>'.( 24) 


Truncating the expansion (1221) at some l = Z max smooths 

small-scale angular fluctuations. _ 

We employ the A-body scheme described by ISellwoodl 
ll2003h that uses a ID grid for the radial coordinate 
while retaining the exact angular expansion for each 
particle. The radial nodes are spaced logarithmically, 
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and subdivide a spherical volume into shells of gradually 
increasing thickness. The functions '$!? m are tabulated 
at shell boundaries and evaluated at other radii by lin¬ 
ear interpolation. To ensure radial continuity, individual 
particle masses are divided between neighboring grid 
points according to the cloud-in-cell procedure, which 
also distinguishes between mass interior and exterior to 
a field point. The particles are smeared in angle because 
the expansion (1221) is limited by the adopted ( max . The 
end result may be pictured as each point particle being 
replaced by a smeared density distribut ion wedg ed betw een 
neighboring spherical radial shells. See ISellwoo3 (120031 ') or 
http: //www .physics . rutgers . edu/~sellwood/manual .pdf 
for a full description. 


3.3 M2M Binning 

A particle within the fcth radial shell, which spans the range 
(ffc-i ,r k ), contributes to ((max + 2)(Z max + 1) real terms 
If we were to use the gravity grid shells to define 
the cells in eq. © as in DL07, then n r +1 radial nodes would 
yield a total of n r ((max + 2)(( max + 1) separate kernel mo¬ 
ments lea. 1281) : note that each particle contributes to 

a single radial shell k and to multiple (l, m , a) components. 

Since n r typically ranges in the hundreds, we super¬ 
impose a coarser grid - the M2M grid - consisting of 
k r + l n r nodes, where all nodes in the M2M grid co¬ 
incide with nodes of the finer gravity grid. Coarser radial 
spacing not only reduces the particle shot noise in each bin 
but also improves the efficiency in solving the FOCE by re¬ 
ducing the number of kernel moments to be computed. 

It may seem that one of these grids is superfluous, for 
if the resolution of the M2M grid is adequate to resolve 
the mass density, then it should be adquate to resolve the 
gravitational forces. However, as argued above, the FOCE 
averages over a large number of particle trajectories making 
coarser M2M binning acceptable, while the finer gravity-grid 
is needed to resolve forces and determine accurate particle 
trajectories. 


3.4 Unequal particle masses 

Even though they are less closely spaced than the gravity 
nodes, the M2M grid spacings should still be small, espe¬ 
cially in the inner half of the halo, in order to represent the 
targeted shape. But with equal mass particles, fine binning 
causes larger statistical fluctuations in the cell occupancy. 
We therefore employ particles that are less massive near 
the center and have gradually increasing masses towards the 
outer edge. To achieve this, we use the dimensionless func¬ 
tion W ( L) = Lo + L to weight the initial particle masses. 
Here, L (= \L\) is the specific angular momentum (in our di¬ 
mensionless units, see H4.3I) and Lo a convenient constant. To 
maintain the desired initial density, we select particles from 
the DF weighted by W~ 1 (L) to ensu re that lighter p articles 
are proportionately more numerous (jSellwoodll200a l. 

We store the initial value Wio = W(Li) for each parti¬ 
cle, which we subsequently adjust during M2M as required 
by the FOCE (eq. 1151) . Furthermore, we adopt the initial 
weights Wio for the priors Wi in the entropy term m- 

Since Mt — f f f d 3 vd 3 x, the fiducial particle mass 


m p = 


1 

N 




(25) 


so that the actual mass of each particle is m p times its cur¬ 
rent weight t Vi(t), and m v ]T\ Wio = Mr- 


3.5 Target Moments and Kernel Computation 

Estimates of the local density from a system of N particles 
are notoriously noisy, so we work with integrated masses. As 
in DL07, we define the (l, m, a) mass harmonic component 
in the fctli bin of the M2M grid as 

Mk;im= f r 2 pTm(r)dr, (26) 

dr k _ i 

where values of pf m (r) are derived from eq. (1241) with T(a:) 
replaced by pr(x). 

Applying eq. (1261) to a discrete particle density yields 
the corresponding harmonic mass components m k . lrn due to 
particles 


( m %lrn\ 


N 

= m p F(Qi) 

i,k 


( cos 

sin m(j>i I ’ 


(27) 


where m v is given by eq. (|25|) . Similarly, the kernel JCk(Xi) 
in eq. © is given by 


(K.%. tlm (Xi A 

\fcMm(*0 ) 


m v np(0i) 
0 


( cos mcj>i 
sin mcf>i 


For the FOCE, the deviation becomes 


A£, m = 


m k-,lr 


~ Mu.,, 


Xi in bin k 


(28) 


otherwise. 


(29) 


§3.9 describes how we estimate the denominator 


3.6 The FOCE 

We use a forward difference to integrate the FOCE, a first 
order ODE. The parameter e determines the rate of change 
of Wi, which is also driven by contributions from all other 
particles in the same radial bin. By iterating multiple uf 
times at fixed t with the adjusted constant c/hf until the 
solution of the FOCE is well converged, we improve the over¬ 
all convergence properties of the algorithm and diminish the 
fluctuations in A*. After each intermediate iteration, we up¬ 
date the right side of the FOCE and renormalize to preserve 
the total mass m p JT Wi(t). 

DL07 suggest the value of e should be a small parame¬ 
ter, eo, times [max^t, |A5fc(i)A*,/cr fc |] — 1 . We adopt their sug¬ 
gestion, but since the maximum value is subject to large 
swings, especially at the outset, we limit its time variability 
by employing a moving average. 


3.7 Halo Initialization 

Our method to prepare the model for M2M integration is 
similar to that described by D09, whereas DL07 employ the 
squeeze method to obtain the initial equilibrium, triaxial 
figure. We proceed as follows: 
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(i) Determine the DF for a spherically symmetric, 
isotropic halo of density p s (r) by Eddington inversion, and 
then redefine the DF to be zero for E > E max , where 
-Emax = ^(rmax), with r max = a, the semi-major axis of 
the target ellipsoid. The density tapers smoothly to zero at 
r m ax, even though the truncation in E is sharp. We apply 
the weighting scheme d escribed in H 3.4 1 and select N m par¬ 
ticles in the usual way (ISellwood fe Debattistall2009l . §2.3). 
This procedure yields a self-consistent spherical model that 
is close to equilibrium; the change in central attraction near 
the outer edge caused by the tapered density creates only a 
very mild imbalance. 

(ii) Apply the scaling transformation, eq. ( 1211 ) . to the spa¬ 
tial coordinates of the particle phase-space generated in step 
(i) in order to deform the sphere onto a self-similar ellipsoid, 
without changing the particle velocities. 

(iii) Evaluate the potential energy Wjj and kinetic energy 
Kij tensor components using the frozen gravitational field 
of the target ellipsoid (determined as described below) and 
velocities from the spherical equilibrium model, respectively, 
and rotate the phase-space coordinates slightly so that the 
kinetic and potential energy tensors become rigorously di¬ 
agonal. We also adjust the velocity components of each par¬ 
ticle so that on average Wn/2Ku = 1 for all three Cartesian 
components. 

(iv) Evolve the system for 25 dynamical times, defined in 
H4.3l below. in the fixed gravitational field of the target ellip¬ 
soid to allow for relaxation during which the model becomes 
slightly less flattened. 


3.8 M2M grid 

In determining the nodes of the M2M grid, we use the spher¬ 
ical radii of the Nm particles generated in step (i) above to 
ensure that each radial bin initially contains at least a preset 
number of particles, or extends over a maximum prescribed 
number of gravity grid nodes, whichever criterion yields the 
fewer nodes. Thus, bins near the center where the density is 
high have the radial extent of only a few gravity grid nodes, 
while bins are as large as we allow in lower density regions. 


3.9 Target Ellipsoid 

To construct a target potential that is as smooth as possible, 
we duplicate steps (i) and (ii) of the above halo initialization 
procedure for a system of Nt particles, with Nt Nm- The 
target potential is that of the larger number of particles, 
which suffers from milder shot noise. 

Moreover, we draw multiple, independent samples of 
Nm particles each from the Nt target population in order 
to estimate the mean target moments, M/* im , and associated 
dispersions, (ea. 1291) . based on the statistical spreads 

of the samples. 


4 AN EXAMPLE: THE EINASTO HALO 
4.1 Mass Density 


Navarro et al.l (1 20041 . l20ld l , iMerritt et al.l (120051 . l2006h , and 
Gao et all (|2008l) show that the logarithimc slope of the 


Label (e y ,£ z ) Axis Ratios Triaxiality 


T a (0.60,0.80) (0.80,0.60) 0.563 

T b (0.70,0.80) (0.71,0.60) 0.766 

P (0.80,0.80) (0.60,0.60) 1.000 

Table 1. Axis ratios and triaxiality parameters of our three 
Einasto halo models. 


spherically-averaged mass density in well-resolved cosmlog- 
ical halos can be fitted as a power law, d In p(r)/d In r = 
— 2 (r/r s ) K , with the value of k dependent on the mass of 
the halo. The scale r s is the radius at which the logarith¬ 
mic slope = —2 . The density that ha s this property is the 
Einasto profile (lEinasto fe Haudlll989l ') 


Pv(r) 



}■ 


(30) 


which rises steeply t owards the center, but remains finite. 
ICardone et al.l (l2005l l give formulas for the mass, potential, 
and force in terms of the incomplete gamma function y(p, x) 
and its complement F(p, x), which is distinguished from the 
usual Gamma function F(a:) by the number of arguments. 
The mass profile is 


3 2 (MM 

K AC \r s J J 

which yields a total mass 


M E (r) = Mo(-J — 7 


Me(oo) = M 0 




(31) 


(32) 


and a gravitational potential 


GMe(oo) f7[M(£)"] , r [M(£)"] 1 

r.r(l) \ r/r s + (*/ 2 )V« j ^ 


We adopt this profile for our spherical mass mode l, fix¬ 
ing k = 0.17, which is the value [Navarro et al.l (120041 ') pre¬ 
ferred for their th ree prototype galax y groups and highest 
resolution haloes dNavarro et al.ll2010l ). 

We set r max = 15, so that the total mass of our models 
is 1.497Mo. (Note Ms(r max ) = 1.960 using eq. dlTTl) : the dif¬ 
ference is due to energy truncation that tapers the density 
smoothly to zero at r max .) We select particles with individ¬ 
ual masses using L 0 = 0.1 (see m which causes ~ 95% of 
all particles to lie within r = 8 , a radius we label as ro. 95 . 

From the discussion in E3J the density of our Einasto 
ellipsoid is 

p e (x) = OyzpE(^), (34) 


which obviously reduces to Pe(x) for e y = e z = 0. 


4.2 Three models 

Three ellipsoidal halo models summarized in Table [l] are 
the focus of this work. All have Einasto index k = 0.17, a 
semi-major axis truncated at r max = 15r s , and are strongly 
flattened (c/a = 0 . 6 ) along the minor axis, but differ in 
their intermediate axes. Although their triaxiality parame¬ 
ters, which are typical of cosmologically formed halos (JS02, 
SFC12), increase from model Ta to P, model Ta is the most 
strongly triaxial. 
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Figure 1. The solid curves show the target mass harmonic profiles up to l = 2 for model P. The dashed curves show the corresponding 
profiles of an ellipsoid of density given by eq. m without an energy cut off. 



Aq 

SmP 

Ta 

T b 

P 

r s (kpc) 

20.3 

20.0 

20.0 

20.0 

20.0 

P 200 (kpc) 

177.3 

232.0 

209.9 

210.1 

210.1 

r max (kpc) 

39.5 

31.3 

33.4 

31.7 

30.7 

Umax (km S _1 ) 

203.2 

216.3 

206.2 

210.1 

215.7 

C200 

8.7 

11.6 

10.5 

10.5 

10.5 

M 200 (x10 11 ) 

12.95 

13.48 

9.97 

10.0 

10.0 


Table 2. Rotation curve ch aracteristics o f ha lo m odels. Aq refers 
to halo model Aq-D-2 in iNavarro et alj fi 2 01 il l. The column 
marked SmP refers to the smooth ellipsoid p e with the eccen¬ 
tricities of the P model 

Ellipsoids have reflection symmetry about their three 
principal planes. We therefore align our coordinate rr-axis 
with the major axis and the z-axis with the polar axis in 
order that only the even terms of the surface harmonic ex¬ 
pansion are non-zero. Oriented in this manner, the m = 0 
components of the 1 /2 is 1 harmonics characterize the ec¬ 
centricity along the minor axis, while even m ^ 0 terms 
describe the density along the intermediate axis. We retain 
even l terms up to £ max = 4 only. 

Fig. [T] illustrates the profiles of the three main contri¬ 
butions - monopole and quadrupole moments - for model P 
and compares the l SC 2, computed from the smooth 

density p e , with those for Adi m (r) for the target ellipsoid. 

4.3 Units 

We employ natural units, for which r s = Mo = 
G = 1, throughout. The unit of velocity is therefore 
Ddyn = ( GM/r s ) 1 ^ 2 and the dynamical time is Td yn = 
{rl/GMf/ 2 = r s /udyn. 

We compare our models with the c osmologically- 
simula ted halo model Aq-D-2 extracted from INavarro et al.l 
(2010, their Tables 1 and 2), which also has the Einasto in¬ 
dex (k = 0.17). Fig. [2]shows rotation curves computed from 
the spherically averaged mass, v ro t (r) = [Moo( r )/ r ] from 
which we deduce the peak circular u ma x and virial W 200 ve¬ 
locities. Setting the length scale to r s = 20 kpc, close to that 
of Aq-D-2, and Td yn = 50 Myr, we have M = 6.74 x 10 11 Mg 
and velocity unit Vd yn — 391 km s -1 . The figure highlights 
the impact of spherically averaging the target density on 
otherwise identical models of different eccentricities. 



Figure 2. Spherically averaged rotation curves V c (r) = 
^GMp 0 (r)/rj ' (in natural units) for each model. The lower 
curve (black) is that of the underlying spherical halo, while in¬ 
creasing the triaxiality parameter yields higher velocity peaks: 
Ta (red), Tb (green), and P (blue). 

4.4 Implementation 

Each of the halo models in Table [T] to be simulated in an 
M2M run consists of Nm = 1.2 x 10 6 particles, while their 
associated target models - from which we derive the target 
potential, moments and dispersions - have Nt = 4.5 X 10 7 
particles. 

With the parameter Lo = 0.1, we find ~ 2.5 x 10 5 parti¬ 
cles inside a radius r = 0.1 (first 16 grid nodes) for the sim¬ 
ulation models. To obtain a specihc M2M binning scheme, 
we sample Nm particles from the the target population then 
identify aggregations of gravity grid nodes with a minimum 
of 2.5 x 10 s particles or a maximum of 10 such nodes, ex¬ 
cept near the center where strict adherence to the criterion 
makes the first few bins too wide: we shorten the first bin to 
5 nodes, which contains between 6000 to 8000 particles; con¬ 
versely, at the boundary we widen the last occupied bin to 
20 to 30 nodes so that it contains > 1000 particles. This par¬ 
ticular scheme typically yields a maximum of 60 M2M radial 
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bins (56 of which are occupied initially) to encompass the 
501 gravity-grid nodes. In cases where a given target value 
M% im vanishes due to symmetry, sampling typically yields 
a nonzero mean with a dispersion that greatly exceeds the 
mean; we omit the corresponding A£. im when solving the 
FOCE to further improve efficiency. 

Step (ii) in 93.71 achieves the ellipsoidal arrangement, 
and (iii) a phase-space near equilibrium. The velocity ad¬ 
justments in the last step yield a mild radial anisotropy, 
P (BT08, see also 95.311 that gradually increases from > 0 
(isotropy) at the origin to a maximum near the boundary as 
the triaxiality parameter increases, i.e., p ~ 0.1 for model 
Ta, while P ~ 0.3 for model P. 

We compute accelerations from the rigid target poten¬ 
tial during the M2M evolution, and advance the motion of 
all particles using leapfrog with a time step At = 0.0025. 
After each time step, we solve the FOCE using a fixed mini¬ 
mum number uf of iterations, gradually increasing uf (from 
rip — 5 initially to a maximum of 12) so that at the end of 
the M2M evolution |A£. im | < 0.1 on average. To achieve this 
threshold we choose fx = 0.5 and Eo = 0.005. Although the 
entropy in the P and Tb models converges sooner, we still 
iterate all models for 200 dynamical times, in order to allow 
particles in the outer halo envelope to complete more than 
one full orbit and to enable most off-grid particles to return. 
Only a negligible number remains permanently outside the 
grid, as reported below (Table [4]l. 

After the end of the M2M evolution, we continue the 
evolution for a further 100 dynamical times with a fully self- 
consistent simulation. This continued evolution enables us 
to check the long term survival of the shape imposed by the 
M2M technique. We also increase Z max to 6, to better cap¬ 
ture the strong flattening along the polar axis, and include 
odd-Z (< /max) contributions to the force determination, as 
further tests of the robustness of the shape. 


5 RESULTS 

Tables [3] and [4] summarize the range of values of key M2M 
quantities for our three models near or at completion of the 
procedure. Since the behavior of the M2M objective, C, is 
essentially stochastic, Table[3]gives the best and worst values 
of the objective to which we add the associated entropy 5, 
cost function C and bin threshold, max |A£. im |. 

The more strongly elongated models have greater order, 
a trend that is reflected in AS'(oo). To support this inter¬ 
pretation, we reran a version of model Ta - model Ta2 - 
for which we omitted steps (ii) and (iii) of the initializa¬ 
tion, i.e., used a spherically symmetric halo, but otherwise 
used the same parameters and procedure during the M2M 
evolution as for the Ta model. The resulting entropy con¬ 
verges to the lower value —0.80; hence, nearly 40% of the 
prospective change in the entropy is captured through the 
direct distortion of the spherical model in our initialization 
process. 

Table []] gives the percentage of off-grid particles and 
of particles that have been driven to a zero weight, No(oo). 
These are but a tiny fraction of Nm ■ The percentage of mass¬ 
less particles that lie within the radii ro.95, -/Vo (ro.95), and 
ro. 95 / 2 , JVo(ro. 95 / 2 ), at the end of the M2M evolution indi¬ 
cates that most massless particles are near the outer edge. 


Ta 

T b 

P 

B W 

B W 

B W 


Objective (£) 

- 0.39 

-0.52 

-0.43 

-0.75 

-0.45 

-0.64 

Chi-square ( C ) 

0.27 

0.53 

0.32 

0.96 

0.31 

0.68 

Entropy (A<S) 

-0.52 

-0.52 

-0.54 

-0.54 

-0.59 

-0.59 

max 1 A? , I 

1 fe;(m 1 

0.08 

0.14 

0.09 

0.16 

0.07 

0.13 

Table 3. Best (B) and worst (W) M2M objective 

near completion 

of algorithm. Other quantities associated with the same instances. 



Ta 

Tb 

P 

Off-grid 

0.002 

0.005 

0.025 

(Vo( 00 ) 

0.033 

0.075 

0.097 

No (ro.95) 

0.012 

0.025 

0.035 

No (ro. 95 /2) 

0.005 

0.009 

0.014 


Table 4. Percent of off-grid and zero-weighted particles at the 
end of the M2M evolution. 

The increase in the numbers of massless particles with the 
triaxiality parameter is strongly correlated with the lower 
(greater) entropy change (ordering), whereas the off-grid 
particle count, which also appears to be correlated with 
triaxiality, is actually driven by the initialization method. 
Thus, there are no off-grid particles for model Ta 2 , sup¬ 
porting the earlier observation of the effect of initialization 
steps (ii) and (iii) as tending to generate highly eccentric 
outer orbits. In this case also No(oo) ~ 0.18% - i.e., there 
are five times as many zero mass particles as in Ta- 

Fig. [3] shows the radial variation of M; m (r), or target 
profiles, and mi m (r), derived from the particles, for model 
P. The curves are drawn to r = 10, which contains > 98% of 
the particles and at two times: the end of both the M2M and 
self-consistent runs; their close agreement in the inner halo 
gives an indication of how well the desired shape is achieved. 

For a more quantitative measure, we give values 
of the cumulative relative error, 5i m = |mj m (r. 95 ) — 

Mlm(r.95)\/\Ml m (j'. 95)1 in Table [5] for the principal (Z, m) 
contributions explicitly included in the M2M algorithm. Not 
surprisingly, the relative errors are smallest for the monopole 
terms, which are numerically the largest. Since all three ha¬ 
los are equally strongly flattened, the smallest errors are for 
the (Z, 0) components, while the errors for the (Z, m > 0) 
components decrease with increasing triaxiality parameter, 
a relationship that is reinforced in the prolate model P. This 
trend can also be understood in terms of the sampling un¬ 
certainty due to the relative contributions from the higher 
harmonics: decreasing the triaxiality parameter at fixed flat¬ 
tening decreases the contribution of the m > 0 harmonics 
and increases the associated sampling uncertainty relative to 
the m = 0 harmonics. Note that for the self-consistent runs, 
the smaller error in the m > 0 components as triaxiality 
parameter increases follows from the same reasoning. 

5.1 Testing Self-Similarity 

iKal /j d 199lfl . JS02 and iDebattista et al.l (120081 1 analyze the 
shapes of their halos directly from the particles by fitting 
an ellipsoid over iteratively refined particle volumes, while 
SFC12 diagonalizes the local inertia tensor. 

We avoid these numerically costly techniques and in- 


© 2014 RAS, MNRAS OOO.ITlITtl 















Self-similar triaxial halos 9 



0 

2 

4 6 

8 

10 

0 

2 

4 6 

8 

10 

0 

2 

4 6 

8 

10 



Radius 





Radius 





Radius 




Figure 3. Target (black), M2M (red) and self-consistent (green) mass harmonic cosine curves of the model P. The black curves are 
barely visible due to the excellent agreement of the red curves. 











Radius 


Radius 


Radius 


Figure 4. Monopole and quadrupole density moments of P (top row), Tb (middle row), and Ta (bottom row) models. Curves displayed 
are at the ends of the M2M (red), self-consistent (green) and perturbed (blue) phases. 
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t a t b p 


(l, m) 

M2M 

SC 

M2M 

SC 

M2M 

SC 

(0,0) 

-3.94 

-2.76 

-3.53 

-2.87 

-3.50 

-3.15 

(2,0) 

-4.43 

-1.95 

-3.11 

-2.39 

-3.08 

-2.72 

(2,2) 

-2.27 

-1.25 

-2.40 

-1.39 

-2.48 

-1.42 

(4,0) 

-1.66 

-1.31 

-1.56 

-1.35 

-1.47 

-1.74 

(4,2) 

-2.46 

-0.61 

-2.57 

-0.89 

-3.01 

-1.02 

(4,4) 

-1.56 

-0.19 

-1.56 

-0.94 

-1.89 

-1.45 


Table 5. Log 10 of cumulative error Si m of the mass harmonics 
at end of both M2M and self-consistent (SC) runs. 


stead use a least-squares fit with Chebyshev polynomials of 
the values of at each node k of the gravity grid, but 

require the value of the fit and its derivative to vanish at the 
origin. Differentiation of the fitted function, Pi^^r), yields 
the corresponding weighted-densities, r 2 pf m (r), shown in 
Fig. [4] These components are combined as in eq. (1221) to 
yield an approximate expansion for the total density 

Uax l 

p(r) « EE'*- nr (0) [ cos m(t> Pim (r ) + 

1=0 m =0 

sinm</>p? m (r)], (35) 

where pf m {r) = r~ 2 dpf m /dr. 

Eq. (1351) forms the basis for testing self-similarity. 
We examine both the radial eccentricity profiles and pla¬ 
nar equidensity contours of our halo models by fitting el- 
li pses to the contou r s wit h the least-squares algorithm of 
I Bender fe MollenhoJ (1 19771 ). a procedure that allows eccen¬ 
tricity, orientation and centering to be estimated simultane¬ 
ously. The left-hand column of Fig. [o]shows the eccentricity 
profile of the primary models, while the middle and right¬ 
most columns show, respectively, the density contours on 
the XY- and ZX-planes. Table [S] gives the eccentricity val¬ 
ues at a fixed set of points along the major axis evaluated 
at the ends of both the M2M and unconstrained evolution. 
Of particular interest are those of the self-consistent runs as 
they correspond to freely evolving, equilibirum halos. Ex¬ 
cept for mild deviations in the outer halo, the tabulated 
values show remarkable uniformity along a given principal 
axis. It is noteworthy that the magnitude of these devia¬ 
tions decreases with increasing eccentricity. Thus, while the 
eccentricity along the minor axis is very stable for all three 
models (< 3% error), along the largest departures (~ 5%) 
on the intermediate axis arise in the model with the least 
eccentricity. The centers of all computed ellipses differ by 
no more than a part in 10 3 , while the direction of the major 
axis fluctuates by less than 2°, declining outwards, with no 
systematic difference caused by the unconstrained evolution. 

In order to compute projected surface densities, we note 
that eq. © can be rewritten as £ 2 = q 2 + z 2 / (1 — s 2 ). Thus 
a projection along the 2 -axis yields a surface density a(q). 
For a projection onto the XT-plane, we integrate along the 
polar axis to obtain a surface density expansion 

‘max 

T,(R,tp) = [E))j (J?) cos mv? + T, s m (R) sin nvp] (36) 

m=0 

such that 


Em(fl) = f] 7Z■ 

l=m 



dz nr 



Pim(r), 


(37) 


where r = (R 2 +z 2 ) 1 ^ 2 . Table[7]gives computed eccentricities 
for all three models for the self-consistent runs. To perform a 
similar projection along a different axis, particle coordinates 
are first rotated so that the axis of projection becomes the 
new polar axis, then apply the above equations. 


5.2 Stability 

If there were any linearly unstable modes present in our 
models, they would be seeded by shot noise. But all three 
models retain close to a constant eccentricity and axial align¬ 
ment in their interior over 100 dynamical times of self- 
consistent evolution, hence any such modes must have neg¬ 
ligible growth rates. 

To test for nonlinear stability, we conduct an additional 
simulation in which we introduce into each of our halo mod¬ 
els a perturber that is turned on adiabatically over a 50 
dynamical times period and turned off over the same time 
period; we refer to these as perturbed runs. We try these 
experiments on the models at the end of the respective self- 
consistent run. Our adpoted perturbation is the l = 2 den¬ 
sity components of an Einasto ellipsoid with the same scale, 
same dimension and same ct-parameter as the host halo, but 
only 10% of the mass and different eccentricities. We use 
a perturber with (e y ,e z ) = (0.55,0.85) for the Ta m odel , 
(0.65 ,085) for Tb, and (0.75,0.85) for P (see lMav fc Binnevl 
1986, for a similar concept), and our experiments terminate 
after 100 dynamical times. In figs. [4] and [5] the corre¬ 
sponding results for perturbed runs are shown in blue. Some 
quantitative differences are evident: model P is no longer 
perfectly prolate and deviations from the target self-similar 
profile are < 10% very near the center. Nevertheless, the fig¬ 
ures remain qualitatively unchanged and may be therefore 
deemed as highly stable. 


5.3 Anisotropy 

We compute the spherically-averaged anisotropy parameter 
/3(r) = 1 — [o-g(r) + a 2 (r)\/2a 2 (r) (BT08) over the range 
0 < r ^ ro. 95 . Here, a 2 (r), crj(r) and a 2 (r) are respectively 
the variances of the polar, azimuthal, and radial velocities. 

Fig. El shows anisotropy profiles for our three principal 
models at the end of their respective self-consistent runs, 
but also includes four other models with lower triaxiality pa¬ 
rameters for comparison. These are a sphere (S), an oblate 
spheroid (O) with e z = 0 . 8 , and two low triaxiality ellipsoids 
with (, e y ,£ z ) = (0.3, 0 . 8 ) (To) and (0.45,0.8) (Td). These 
additional models were created in the same manner as our 
principal models through the self-consistent runs, but for 
the sphere we fit only the monopole component during the 
M2M run and admit l m ax = 2 for the self-consistent evolu¬ 
tion. Fig. [ 6 ] confirms that the sphere has /3 « 0 as expected, 
especially in the well-resolved interior region. The tangential 
bias (/3 < 0 ) near the boundary is an artifact of enforcing the 
targeted monopole moments in that region; the same effect 
can be seen in the oblate and triaxial models Tc and Td- 
Flattening causes a radial bias (/? > 0) in the inner regions, 
which is enhanced in the mildly triaxial models, gradually 
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Figure 5. The left panels show the eccentricity profiles derived from least-squares ellipse fits to P (top), Tb (middle), and Ta (bottom) 
models. The middle and right panels of each row show density contours on the XY - and XZ- planes, respectively. The colours in all 
panels refer to the ends of of the M2M (red), self-consistent (green) and perturbed (blue) evolution. 


overcoming the tangential bias at the boundary. Indeed, the 
fiducial models, which have the highest triaxialities, are ra¬ 
dially biased throughout. 

The thick black lines in the sa me fi gur e are a dapte d 
from the middle panel of Fig. 3 of Ludlow et al.1 (1201 ill . 
which consists of halos in the Millennium-II simulation with 
Einasto parameter a = 0.178. The broken black line shows 
the median anisotropy profile of their halo sample, while the 
thick black lines represent the upper and lower bounds of an 
envelope that encloses the bulk of the plotted profiles. Pro¬ 
files for our fiducial models all lie within the envelope, while 
those of the added models with low triaxiality either straddle 
or lie beneath the lower boundary, a result consistent with 
the high triaxiality, prolate bias exhibited by cosmological 
halos (SFC12). Thus, even as an idealization, the self-similar 


triaxial geometry is supported by an anisotropy profile in¬ 
distinguishable from that of cosmological halos with compa¬ 
rable triaxiality, a noteworthy result given that we have not 
constrained the velocities in our models. 

5.4 Run time 

In keeping with the averaging property of the FOCE (see 
& every particle ought to complete several orbits to en¬ 
sure convergence. Periods of orbits with radius near and 
beyond r 2 oo of Einasto halos easily exceed 100 dynamical 
times so strict adherence to such a criterion would require 
excessively long running times. However, a more reasonable 
length of M2M evolution achieves robust self-similarity in 
the halo’s interior while self-similarity in the outer region 
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t a t b p 



= 

0.60 

= 

0.80 

e y = 

0.70 

£z = 

0.80 

£ y = 

0.80 

£z = 

0.80 

X 

M2M 

SC 

M2M 

SC 

M2M 

SC 

M2M 

SC 

M2M 

SC 

M2M 

SC 

0.10 

0.585 

0.613 

0.816 

0.809 

0.703 

0.700 

0.798 

0.791 

0.811 

0.808 

0.816 

0.806 

0.15 

0.581 

0.604 

0.808 

0.798 

0.698 

0.694 

0.791 

0.784 

0.803 

0.798 

0.807 

0.800 

0.25 

0.583 

0.596 

0.799 

0.786 

0.695 

0.689 

0.791 

0.779 

0.796 

0.788 

0.799 

0.795 

0.50 

0.599 

0.598 

0.796 

0.793 

0.702 

0.698 

0.803 

0.790 

0.796 

0.789 

0.798 

0.804 

0.75 

0.603 

0.602 

0.800 

0.801 

0.703 

0.704 

0.805 

0.800 

0.800 

0.796 

0.801 

0.812 

1.00 

0.599 

0.607 

0.801 

0.803 

0.700 

0.702 

0.802 

0.802 

0.799 

0.795 

0.800 

0.808 

2.00 

0.606 

0.613 

0.803 

0.814 

0.696 

0.698 

0.800 

0.798 

0.799 

0.797 

0.799 

0.792 

3.00 

0.599 

0.602 

0.790 

0.794 

0.698 

0.698 

0.799 

0.794 

0.798 

0.796 

0.794 

0.796 

4.00 

0.608 

0.629 

0.797 

0.807 

0.701 

0.714 

0.798 

0.800 

0.799 

0.798 

0.796 

0.801 

5.00 

0.599 

0.621 

0.800 

0.807 

0.701 

0.700 

0.801 

0.786 

0.798 

0.793 

0.798 

0.797 

6.00 

0.597 

0.618 

0.786 

0.804 

0.693 

0.691 

0.799 

0.787 

0.799 

0.790 

0.799 

0.792 


Table 6. Eccentricity profiles at ends of the M2M and self-consistent (SC) evolution. 


Ta T b P 


X 

Ey = 0.60 

Ey = 0.70 

Ey = 0.80 

0.10 

0.618 

0.678 

0.807 

0.15 

0.601 

0.674 

0.795 

0.25 

0.582 

0.679 

0.783 

0.50 

0.567 

0.701 

0.788 

0.75 

0.572 

0.707 

0.794 

1.00 

0.591 

0.703 

0.793 

2.00 

0.611 

0.699 

0.787 

3.00 

0.600 

0.696 

0.786 

4.00 

0.612 

0.688 

0.796 

5.00 

0.611 

0.680 

0.790 

6.00 

0.614 

0.687 

0.790 


Table 7. Profiles of the eccentricity of the surface density pro¬ 
jected onto the XY-plane at end of the self-consistent runs. 

of the halo remains less than perfect. This is illustrated in 
Fig.Qwhich displays the quadrupole densities p 2 o{r) (right) 
and p 22 {r) (center), and the eccentricity profile (left) after 
running times of 100 (red), 150 (green), 200 (blue) and 250 
(cyan) dynamical times and the subsequent self-consistent 
run for the Ta model. We have omitted the monopole den¬ 
sity poo since all three cases align perfectly with the targeted 
density, which suggests it is the first moment to fully con¬ 
verge. It is evident that the P20 moment converges sooner 
than p 22 , especially near the boundary, which is reflected 
in the eccentricity profile. We have verified that model P 
has converged adequately by 150 dynamical times. Table 0 
underscores the connection between accuracy of fitting the 
various moments, high eccentricity and its reduction in the 
uncertainty of the moment due to particle noise. In our im¬ 
plementation, more eccentric models converge more rapidly. 


6 SUMMARY AND CONCLUSIONS 

We have used the M2M method to create collisionless halo 
models with self-similar, ellipsoidal geometry and a mass 
density profile relevant to cosmological halos. Our model 
halos are self-similar, which is not in strict accord with cos¬ 
mological results (JS02, SFC12), but is a reasonable approxi¬ 
mation with theoretically advantageous properties. They are 
the first step in an ongoing study of the collisionless growth 
of stellar disks within triaxial halos. 



Figure 6. Anisotropy profiles arranged in ascending triaxiality 
parameter: Spherical (thin black line), Oblate (red), To (green), 
Tq (blue), Tb (cyan), Ta (magenta), and P (yellow) halo models 
(see text for definition of non-fiducial models). Superimposed on 
these are the upper and lower b oundaries (thick black lines) of the 
profiles in Fig. 3 of iLudlow et al] d201 ill and the corresponding 
median curve(broken black line). 

Both DL07 and D09 have previously created triaxial 
halos by this method, but there are minor differences in our 
approach, and we have adopted a more cosmologically moti¬ 
vated model. The main difference in methodology is that we 
do not implement a moving average, (eq. HGII . which we find 
actually degrades the results and we show why it is superflu¬ 
ous. Other minor differences in our approach are described 
in ilTTl 

DL07 introduced the M2M objective as a goodness-of- 
fit measure to accommodate observational data. Our present 
application is to a theoretical target model, and we high¬ 
light the similarities and differences between M2M and 
Schwarzschild’s approach. We also follow DL07 and D09 
by interpreting the dispersions as sampling error consistent 
with shot noise inherent in any A-body representation, but 
construct our target potential from a static A-body model 
with a very large number of particles, which we use not only 
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Figure 7. Quadrupole weighted densities in = 0 (left), m = 2 (center) and eccentricity profile of the Ta model for target (black) and 
M2M run times of 100 (red), 150 (green), 200 (blue) and 250 dynamical times (cyan) shown after respective self-consistent runs. 


to compute the target potential, but also to determine the 
mean and dispersions that serve as proxies for Pj and Uj, 
respectively, in eq. 

The full procedure consists of several distinct steps: af¬ 
ter initializing the particles as described in we allow 

the model to relax in the smooth target potential for 25 dy¬ 
namical times. We then apply the FOCE over 200 dynamical 
times, a sufficient time span to ensure that changes in the 
entropy level are minimal. 

Our adpoted model is the Einasto halo and we present 
three aspherical models that bear strong simil arities to halos 
from the Millenium-II simulation selected bv lLudlow et al.i 
ll201ll) . We demostrate that they are stable by evolving them 
in their self-consistent potential over a period of 100 dynam¬ 
ical times, during which their properties barely changed. We 
also tested non-linear stability by perturbing the model and 
showed that its shape was restored quite well after the per¬ 
turbation was removed. In additional experiments, not pre¬ 
sented here, we have found that singular mass profiles, such 
as the NFW or Hernquist, present no difficulty, and higher 
eccentricities can also be achieved, with the proviso that 
Imax may need to be increased, which would lengthen com¬ 
putation times. 
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